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ABSTRACT 

Acceleration methods are presented for solving the steady state 
incompressible equations* These systems are preconditioned by introducing 
artificial time derivatives which allow for a faster convergence to the steady 
state* We also consider the compressible equations in conservation form with 
slow flow* Two arbitrary functions a and 3 are introduced in the general 
preconditioning. An analysis of this system Is presented and an optimal value 
for 6 is determined given a constant a. It is further shown that the 
resultant incompressible equations form a symmetric hyperbolic system and so 
are well-posed. Several generalizations to the compressible equations are 
presented which generalize previous results. 
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1 • INTRODUCTION 


In this study we consider ways of reaching a steady state for the 
incompressible fluid dynamics equations and also for low Mach number 
compressible flows. We shall only consider time-marching schemes that are 
represented by hyperbolic systems, Chorin [6] developed the artificial 
compressibility method which is further discussed by Peyret and Taylor [10], 
We consider generalizations of this method by allowing artificial time 
derivatives in all the equations and not just the continuity equation. This 
allows for faster convergence and also facilitates a uniform treatment for 
both primitive variables and conservative variables. It is shown that the 
resultant equations forms a symmetric hyperbolic system and so is well-posed 
for both primitive and conservative formulations. 

We next consider compressible flow with very low Mach numbers. As is 
well-known, this system is stiff due to the large ratio of the acoustic and 
convective time scales. A number of people have considered preconditionings 
of these equations in various special cases, e.g., Viviand [20], Briley, 
McDonald, and Shamroth [2], Choi and Merkle [4], [5], Rizzt [13], and Turkel 
[17], [19]. In this study we generalize these various approaches. In all 
cases we consider primitive variables p, u, v plus an additional variable. 
After the analysis is complete It is shown how one can reformulate the system 
in conservation form. 

As pointed out by Briley et al., [2], it is necessary to nondiraensionalize 
the equations so that the pressure does not go to infinity as the Mach number 

A 

goes to zero. In [2] this is accomplished by choosing (p ,Cp,u,v,hg) as the 
dependent variables. In this study we shall concentrate on (p,u,v) 
variables or else the conservation variables (p,pu,pv,E). Instead in the 
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Runge-Kutta code [8] used the variables are nondiraensionalized so that p = p 
= 1 in the far field. 

2. INCOMPRESSIBLE FLOW 

In this section we consider the incompressible inviscid equations. 
Extensions to the viscous equations will be considered later while the 

following sections will discuss the effects of compressibility. We will only 
consider time-independent solutions. Nevertheless, since we shall discuss 
time-marching algorithms we begin with the time-dependent incompressible 
inviscid equations • 

u + v =0 

x y 

u + uu + vu + p = 0 (2.1) 

t x y r x 7 

v. + uv + vv + p =0. 

t x y y 

These equations can also be written in conservation form as 

u + v =0 
x y 

u t + (u 2 + p) x + (uv) y = 0 (2.2) 

v ^ + (uv) x + (v 2 + p) y = 0. 
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In this study we shall only consider smooth solutions to the systems (2.1) and 
(2.2). The only discontinuous solutions of interest are contact 
discontinuities, vortex sheets, etc., which are essentially linear phenomena 
and so are extensions of smooth flows. Shocked flows are not of interest for 
these equations and hence, the systems (2.1) and (2.2) are identical. 

Since we are only interested in steady solutions we will modify the time 
derivatives that appear in (2.1) and (2.2). The simplest such modification is 
the pseudo-compressibility approach which adds a pressure time derivative to 
the continuity equation [6], [10], [15]. Then all the equations can be 
marched in time until a steady state is reached. We shall consider 
generalizations of this technique. All the time-dependent equations that we 
consider form hyperblic systems. Since there is no decay mechanism except for 
boundaries we can accelerate to a steady state only by increasing the 
allowable time step. By normalizing the fastest speed it is shown in [19] 
that we accelerate the convergence when all the speeds are close together in 
absolute value. Conversely, the worst convergence occurs when the speeds of 
the system differ by orders of magnitude. It is also shown in [19] that in 
order to have a well-posed problem that is compatible with the steady state, 
especially in terras of boundary conditions, it is desirable to have a 
symmetric hyperbolic system. For a symmetric hyperbolic system, when the 
preconditioning matrix is positive definite we are guaranteed that we have not 
changed the appropriate number of boundary conditions and that we have not 
introduced any nonphysical time reversals. 
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We therefore consider the following extension of system (2.1) 


— p + u + v 
0 2 t x y 


= 0 


au , . 

— r- + U + UU + VU + P 

0 2 t t x y *x 
p 


= 0 


(2.3) 


av 



+ v. 


+ uv 4- w + p 
x y *y 


0 . 


Here, a and 3 are functions to be determined. When, a = 0 we recover the 
standard pseudo-compressibility method and we need only determine 3. To form 
a conservation system we multiply the first equation by u and also v and 
then add to the second and third equations respectively. The resulting system 
is 


— p +u + v = 0 
B 2 t x y 


■ ( ” -- y )U P t + U t + (u2 + p) x + (uv) y = 0 


(2.4) 


— 2 ~ " P t + V t + ^ UV ^x + + P V = °* 


B 


Note 1 : The system (2.4) is not truly conservative for time-dependent 
flows. However, we have in any case destroyed the time accuracy and the 
system is fully conservative in the steady state. 

Note 2 : Even for the original pseudo-compressiblity approach, a = 0, one 
should add pressure time derivatives to the momentum equations in the 
conservation form. Some authors, e.g., [4], [13] have not added these 
derivatives which amounts to choosing a = -1. 
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Note 3 : We shall do all the analysis on the nonconservative system 

(2.3) . All the results will be equally valid for the conservative system 

(2.4) with the appropriate a and B. 

We first rewrite (2.3) in vector form as 



or 


with 



w = (p.u.v) 1 . 


Multiplying (2.6) by E we rewrite (2.5) as 



( 2 . 8 ) 
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with 

A — EA Q y B = EBq • 

In order to consider the wave speeds of (2.7) we Fourier transform the 
system. The wave speeds of (2.5) are given by the eigenvalues of 

D = Uj A + u 2 B, -1 < Up u 2 £ 1, (2.9) 

where iiij, are the x and y components of the Fourier transform 

variable. Defining 

q = uu ^ + vu 2> (2.10) 

we find that the eigenvalues of D are 


and 


d 0 “ 9 , 


(2.11a) 


d . = l h 


(1 - a)q ± )/o. - a)^ q^ + 43^ 


(2.11b) 


Remark : In the special case a = 1, we have = ±3 and so the 

"acoustic" sound speed is isotropic and independent of the flow velocity. 

We note that for all values of a and 6, d + and d_ have opposite 

2 

signs, i.e., d + • d_ = -6 is always negative. This corresponds to subsonic 
flow for a compressible fluid which is appropriate for the incompressible case 
being considered. 

We next consider the choice of 3. We consider a as given and we wish 
to choose 3 to minimize the largest possible ratio of wave speeds. Thus, we 
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wish to choose 3 so as to minimize max( | d^/dj | ) i,j = 0, ±. After some 
algebra we find that the appropriate 3 is given by 



a < 1, with condition no. |2 - a| 
a > 1, with condition no. a . 


( 2 . 12 ) 


Formula (2.12) is not useful since q, given by (2.10), is a function of the 
Foruier variable (w^,^) while 3 must be given in physical space. Hence, 
we replace (2.12) by 

1 2 - a, a < 1 (2.13a) 

a, a >_ 1. (2.13b) 

The ratio of the fastest to the slowest speed now also depends on the ratio 
2 2 2 

(u + v )/q and will be larger than that given in (2.12) unless 

2 2 2 
q = u + v . 

Remark : It follows from (2.12) that the optimal a is a = 1 in which 

case the condition number is one, i.e., all the speeds have the same 
magnitude. Since 3 cannot be a function of the Fourier variables we must 
use (2.13) which means that the condition number is a function of 
Nevertheless, this is still the best result for a range of Fourier modes In 
multidimensions. We remind the reader that the original artificial 
compressibility corresponds to a = 0 for the primitive variables and to 
a = -1 for the conservative variables. 
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2 

We note that in all these formulae 3 is not constant but rather is a 

2 2 

function of the speed, u + v • To avoid difficulties near stagnation points 
(2,13) must be modified so that 6 cannot approach zero. For example, (2.13) 
can be changed to 


max[ (2 - a)(u 2 + v 2 ),e] 

a < 1 

K max[a(u 2 + v 2 ),e] 

a > 1 


(2.14) 


On dimensional grounds e should be a fraction of (u 2 + v 2 ) . From later 

max 

considerations, (2.16), K should be chosen slightly larger than one. 

Until now we have only discussed the wave speeds, i.e., the eigenvalues 
of D, (2.9). We have shown that these eigenvalues are always real and so 
(2.5) is a hyperbolic system. We next wish to find out whether the system can 
be symmetrized. Gottlieb and Gustaffson [7] have suggested a general 
technique to check if a system can be simultaneously symmetrized. A necessary 
condition is that A and B can each be separately symmetrized. Since A 
can be symmetrized it can also be diagonalized. Furthermore, the 
diagonalization of A is unique except for exchanges of rows and columns and 
also an additional similarity transform using diagonal transformations. Thus 
after A has been diagonalized one only need check if B can be symmetrized 
by a diagonal similarity transform. 

We first need the eigenvectors of A. This is also useful for 
constructing characteristic boundary conditions. If follows from (2.11) that 


a Q = u 


a = 

± 


_ (1 - a)u ± >/( 1 - a ) ^ u 2 + 43 ^ 
2 


(2.15) 
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Then det(T ) = /(I - a)^ u + 4f5^ * 0 and so the transformation is 

nonsingular. Furthermore, the columns of T -1 are the eigenvectors of A. 


It then follows that 



Let Dg = diag(dj ,d^ ,dg) with 

d 1 = a_/(a + - a_)(u - a_), 



d 2 = a +^ a + ” a _K a + “ u )> 
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d 3 = (a + - u) (a_ - u) 


" a + a - 


2 2 2 
3 - a(u 4* v ) 


Then D Q TAT -1 d" 1 is still diagonal, while D Q TBT _1 D” 1 is now 

symmetric. It follows from the definition of a + , a_ that a + > 0 > a_. It 

can then be shown that Dq Is real and hence the system Is symmetrizable if 
2 2 

and only if q _< d + for all (to , 0 ^ 2 ) or equivalently 


q 2 v f 2 , 2 v 

3 > ct(u + v ). 


(2.16) 


We note that from (2.13) we have, that the optimal 3, for a 1, is gotten 

by choosing an equality in (2.16) rather than an inequality. Hence, if one 

wishes the system to be both close to optimal and symmetrizable we should 
2 2 2 

choose 3 slightly larger than u + v . Furthermore, for a < 1, (2.13a) 
implies (2.16) automatically. For a 0 (2.16) is always satisfied for all 

3. 

When using an explicit method we need an upper bound on the largest eigen- 
value of D. A typical explicit scheme has a stability criterion of the form 



(2.17) 


where A is a typical mesh length and K is a constant that depends on the 
scheme. Using (2.11b) we replace d + in (2.17) by Its upper bound 


d + < 


, 2 . 2, 

(u + V ) 


(l-a+V(l-a) 


2 2 
u + v 
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2 2 2 

with 8 / (u + v ) given by (2.14). If we use a general curvilinear mesh 
then the corresponding formula is given by (3.20), with c = °°. 

The previous discussion has been scheme independent and relied only on 
equalizing the wave speeds for the differential equation. We now discuss the 
implementation for some difference schemes. For an explicit scheme the time 
step is restricted by the fastest moving wave. Thus, the previous analysis 
insists that the time step chosen by a stability analysis should not be 
inappropriate for the slower waves. If the wave speeds differ significantly 
then the slower waves will propagate very slowly and convergence will also be 
slow. Furthermore, for most explicit schemes the damping of the scheme is 
small for small At and so the slowly moving waves will not be damped very 
much. Hence, our analysis is certainly appropriate for standard explicit 
schemes. 

Using an implicit method it is less clear that the stiffness of the system 
matters. If one uses a backward Euler method then one can show [8] that for 
large At that one approaches the classical Newton-Raphson iteration 
scheme. In this case the convergence is not very affected by the stiffness of 
the system. In [4] computations are presented that show fast convergence for 
one-dimensional problems. However, in multidimensions It Is not practical to 
invert the matrix that one gets using a fully implicit scheme. Instead one 
frequently uses an ADI type algorithm. In this case one should not choose 

very large At [16] but rather one close to the explicit Courant condition. 

2 

This occurs because of the (At) AB term that is created by the splitting. 
Hence, again a At that is appropriate for the fast waves is inappropriate 
for the slow waves. Hence, our preconditioning which will equilibrate the 
wave speeds will also accelerate ADI type methods. Using the notation of Beam 
and Warming [1] we write an Implicit scheme for (2.4) as 
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We rewrite (2.17) as 


E "‘ l 1 + 4t ( E " h A S + E " h B S)l 4 “ ' ( f x + g y) • <2 - 18) 

We now apply an approximate factorization to (2.18) and ignore errors in the 
conservation form of the left-hand side to get 

Ef 1 ! 1 + At EA 0 ] [i + At Ip; EB q ] Aw - -At(f“ + g^). C2.19) 
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Since the matrices E, A = EAq and B = EBq are well conditioned there is no 
way that the splitting error can slow down the convergence compared with the 
standard ADI splitting. 

For a = 0 we need only invert 2x2 blocks. For general a we can 
use the factorization suggested by Pulliam and Steger [11]. Hence, 



can be diagonalized, i.e., UAU * = D^, VAV * = D^. Ignoring, again, 
conservation errors in the left-hand side of (2.19) we rewrite (2.19) as 


E 1 U 


I + At 


8x °1 


if 1 V 


I + At 




V * Aw = -At(f n + g n 

\ x °y 


( 2 . 20 ) 


and so we need only invert scalar tridiagonal matrices rather than block 
tridiagonal matrices. 

In practice one usually solves the viscous equations rather than the 
inviscid equations. The easiest remedy is simply to add the viscous terras to 
the right-hand side of (2.19). One usually finds, for large Reynolds numbers, 
that the time step is restricted only by the inviscid terms. Hence, there is 
no need to include a viscous Jacobian on the left-hand side of (2.19). 
Furthermore, the preconditioner, E, still equilibrates the inviscid time steps 
and reduces the splitting error In (2.19). 

We next consider the implementation of the scheme on a staggered mesh. 
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V 



V 

Figure 1 


The steady state equations are independent of At and so we retain the 
improved accuracy of the staggered grid independent of our treatment of the 
time-marching algorithm. Thus, for example, we discretize the x-moraentum 
equation in (2.3) by 



and is a function of a given in (2.14). Using an explicit scheme p t 
at (i,j) is already known from the first equation. With an implicit scheme 
we now have contributions of p t in the momentum equations which contribute 
to both the diagonal and off diagonal blocks. 
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In this section we have only considered Cartesian coordinates. The 

extension to general curvilinear coordinates is straightforward. This will be 
done in the sections on compressible flow which will contain the 

incompressible case as a limiting solution. Here we shall only show that the 
matrices are simultaneously symmetrizable in curvilinear coordinates whenever 
they are symmetrizable in Cartesian coordinates. Consider the equation in 
Cartesian coordinates (X,Y), 

w^ + Aw x + BWy =0. (2.22) 

Let x = x(X,Y), y = y(X,Y) be general coordinates then 

w^ + A, w + B, w =0 (2.23) 

t lx 1 y 

with 

A x = Ax x + Bx y , B x = Ay x + By y . (2.24) 

Since A^ and Bj are linear combination of A and B it follows that 

whenever A and B can be symmetrized simultaneously so can A^ and B^. 


3. COMPRESSIBLE (p,u,v,S) SYSTEM 

In the previous section we have considered incompressible flow where the 
unknowns are (p,u,v). We next consider the compressible equations 
concentrating on low speed flow. Since our analysis is local we need only 
consider flows that locally have a small Mach number. The flow can even be 
supersonic in other regions. Hence, it is useful to consider the conservation 
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forra of the equations. In considering the compressible equations we need an 
additional unknown. Three possibilities are entropy S, or density p or 
else to use Bernoulli s law stating that the total enthalpy is constant, i.e., 
isoenergetic flow. In all cases we shall ultimately cast the equations in 
conservation form but the three possibilities lead to different 
preconditioning. As before we shall do the analysis on the primitive 
variables and only at the end shall we derive the conservation variable 
version. In this section we consider the (p,u,v,S) system while the other 
possibilities are discussed in later sections. 

The time-dependent Euler equations can be written in Cartesian 
coordinates (X,Y) as 


~^2 P t + “T (up X + Vp Y } + U X + V Y = 0 

U t + UU X + VU Y + p x /p = 0 

V t + UV X + VV Y + P Y^ P = ° 


( 3 . 1 ) 


S t + uS x + vS Y = 0 

where 

p = p (p , s) . 


We now introduce curvilinear coordinates x = x(X,Y), y = y(X,Y). The Euler 
equations in (x,y) coordinates are 
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— — 7 T p_ + (Up +Vp)+u Y - v X - u Y + v X = 0 

2 r t 2 *x r y xy xy yx yx 

pc pc J J J J J 


Ju + Uu + Vv + (p Y - p Y )/p = 0 
t x y x y *y x 


Jv + Uv +Vv +(-p X +p X )/p = 0 
t x y x y *y x 


JS„ + US + VS =0 
t x y 


(3.2) 


where p = p(p,S), and 


U = uY - vX , V = -uY + vX , J = X Y - X Y . (3.2a) 

y y ’ xx’ xyyx 


We precondition this system by a generalization of (2.5). We thus obtain 
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Since the entropy decouples on the matrix level (we freeze coefficients and so 
p no longer depends on S) for stability theory we can reduce (3.3) to the 
simpler equivalent form 



We note that (3.4) is very similar to (2.5). In fact, setting c = 8 and 
p = 1 and using Cartesian coordinates, (3.4) reduces to (2.5). As before we 
rewrite (3.4) as 

JE _1 w + A« w + B n w =0. (3.5) 

t 0 x 0 y 

Multiplying (3.4) by E we obtain 



0 


(3.6) 
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A = EA 0 = 


3 2 U 


-auU 


pc 


2 p 


-avU 


pc 


P6 2 Y y 


-auY + U 

y 


-avY 

y 



? 

B V 


2 

-pB Y 


B = EB q = 


-auV x 

2 ~ p 

pc 


auY + V 
x 


- otvV + "x 


— avY 


pc 


x 


2 

pr x 


-ctuX 


x 


-avX + Vi 
x 


(3.7) 


and the Jacobian J is given by (3.2a). To find the wave speeds we again 
examine the eigenvalues of 


We define 


and 


D = A + “1 .< w 2 — 


= Y w, - Y , 
1 y 1 x 2 * 


^2 = " X y “l + X x W 2 


q = Uw^ + Va^ = 


(3.8) 


(3.9) 


• U and V were defined in (3.2a). Then the eigenvalues of D are 


d 


0 


= q 


(3.10a) 


and 
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(3.10b) 


We note that without preconditioning we have 


d ± - q ± / £ 2 + £ 2 c . 

Remark ; If we consider the special case o = 1 + g 2 /c 2 = 1, then 

“ ±3/^^ + &2 “ 1 /c - + £ 2 for low speed flow. Hence, the acoustic 

wave moves with a speed independent of the velocities u, v and this wave is 

isotropic except for grid effects. Also, we note that (3.10) is independent 
of c. 

We also know that 


2 2 
1 *2 


- tXy + + (** + ^ - 2(X x X + V x Y )«, „ 2 < L : 


,2v 2 


(3.11) 


x v + Y l + X ! + Y l + 2 I X x + Y Y I 
x x y y 1 x y xy 1 


For subsonic flow d + and d_ have the opposite signs. In fact 


d + d_ = -( £ 2 + £ 2 - q 2 /c 2 )g 2 < 0 


2 2 2 

whenever u + v < c . For an orthogonal mesh 
expression for L 2 , (3.11) simplifies. 


X X + Y Y = 0 
x y x y 


and so the 
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We next consider the choice for We wish to choose 3 so as to 

minimize the largest possible ratio of the wave speeds, i.e., to minimize the 

maximum of the ratio of the d^s in (3.10). In order to simplify the algebra 

2 2 

we assume that q 2 /c 2 1, i.e., slow flow, we also assume that 3 /c ^ 1» 

The optimal 3 is then computed as (cf. (2.12)), 

a < 1 

+ 0(q 2 /c 2 ) . (3.12) 

a > 1 

2 2 2 

The condition number is the same as (2.12) to within 0((3 + 3 )/° )• 

Similar to the incompressible case, (3.12) is not of immediate use since 

and q all depend on the Fourier variable (oj^,^)* see (3.9). Instead 
we suggest using 

B 2 L 2 I' 

u2 + y2 u 


a < 1 


a > 1 


(3.13) 



where L is given in (3.11). 

We next rewrite (3.3) in terms of the conservation variables (p,pu,pv,E) 

o n 

with E = — ]] j + ~ (u + v ). We then obtain 



+ F x + Gy = 0 


(3.14) 
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where F and G are the standard Euler fluxes in curvilinear coordinates, 
J is the Jacobian, and 


z i = 

8 c 


- (y - +2<3L^» 


c 

2 


z, = 


- 1) f-r-f + ( 

T r e 


2a + 1 1 

2 “ 2^ 2 J 


(3.15) 


= z h + (y - l)a(u^ + v^)/e 4 , h = — 5—r- + - u - 2 t y2 
1 Y - 1 2 


2 . 2 W „2 


Thus, as expected, we recover the correct steady state equations. We can also 
eliminate in (3.14) and obtain equations only in terms of p^, (pu)^, 

(pv) t and Ej_. We then have 


pu 

J(I + Q)| I + F + G 
1 pv | x y 


= 0 


(3.16) 


where I is the identity matrix and 


R Z l -UZj -VZj Zj 


Q = 


p 2 2 

R z 2 u -u z 2 


R z 2 v ~uvz 2 


-uvz. 


z 2 u 


-V z 2 z 2 V 


(3.17) 


R z 3 -uz 3 -vz 3 z 3 
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and R 2 = (u 2 + v 2 )/2. When a = 0 and Q 2 = u 2 + v 2 this reduces to the 
preconditioner found in [18] by a different technique. We note that for a = 0 
the optimal 8 given by (3.13) is g 2 = 2(u 2 + v 2 )/L 2 . Furthermore, we can 
invert I + Q simply to get 

(I + Q) -1 = I - Q. (3.18) 

c 

Because of the structure of Q we can multiply Q times a vector using seven 
multiplications . 

As before, the stability criterion for a typical explicit scheme for 
(3.16) has the form 

^-<K/d + (3.19) 

where K is a constant that depends on the scheme. It follows from (3.9)- 
(3.11) that 




where L is defined in (3.11), is a sufficient condition for stability. For 

2 2 2 

slow speed flow we can ignore all terras of the order (U + V )/c and 

8 2 /c 2 . Also since B 2 = 0(u 2 + v 2 ) by (3.13), we see that At is 

independent of c and depends only on the local velocity. As pointed out 

2 2 

previously the special choice a = 1 + 3 /c simplifies the formulas. We 
then find that 

for a = 1 + 3 2 /c 2 . (3.21) 

c 
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As long as the flow is subsonic the square root is meaningful. 

As in the incompressible case we find that the matrices A and B can be 
simultaneously symmetrized when 

3 2 > cx(u 2 + v 2 ). (3.22) 

In forming the preconditioned system (3.16) we eliminated the pressure 
term p t from (3.14). Since we are not interested in the time-dependent 
solution we can instead eliminate p t> 

u 2 + v 2 P 

(“ — 2 ~~) P t = y - J l + u(pu) t + v(pv) t - E fc . (3.23) 

As before, we need to do something special in the neighborhood of stagnation 
points. This system now solves for (p,pu,pv,E) and so is more similar to 
the incompressible limit. 


4. COMPRESSIBLE (p,u,v,p) SYSTEM 

In the previous section we appended the entropy equation to the 
incompressible (p,u,v) equations and did not precondition the S 
equation. This had the benefit that the entropy equation decoupled and so 
even in the compressible case we needed to only consider three equations (see 
( 3 . 2)— (3 .4) ) . Choi and Merkle [4], [5] have discussed a (p,u,v,p) 

formulation which we now analyze in more detail. 

Again considering curvilinear coordinates x = x(X,Y), y = y(X,Y) the 
Euler equations are (compare with (3.2)). 
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where 


U = uY - vX , 

y y 


V = -uY + vX , 
X x’ 


J = X Y 
x y 


X 

y 


Y . 

x 


We precondition this system similar to (3.3) where again the last equation is 
not preconditioning. Thus, we are now not changing the p equation rather 
than the S equation of (3.3). We then obtain 



(4.2) 
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In order to facilitate comparisons with (3.3) we change variables in (4.2) to 

a (p»u,v,S) system. Formally, we define S = Jtn[p/p Y ] and so 

dp = — (dp - pdS). Substituting into (4.2) we get 
c 
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Comparing ( 3 . 3 ) with ( 4 . 3 ) we see that using the (p,u,v,p) formulation has 
introduced an additional p t preconditioning into the entropy equation. 
Hence, the entropy equation no longer decouples from the previous three 
equations. This complicates the analysis. The advantage of the (p,u,v,p) 
is that it simplifies the preconditioning in conservative variables as will be 
seen later. Solving for (p,u,v,S) t we find 
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or 


Jv/ + Aw 

t X 


+ Bw . 

y 


(4.5) 


Comparing the matrices A and B of (3.6)-(3.7) with that of (4.4)-(4.5) 
the eigenvalues of any linear combination of A and B are unchanged since 
the last column has zeros except for the corner element. Hence, all the 
preconditionings that were considered before in Section 3 are equally 
efficient for the (p,u,v,p) system. In particular, an optimal 3 is given 
by (3.13) and the time step restriction for a typical explicit method is given 
by (3.20) and (3.21). 

We now rewrite our preconditioned (p,u,v,p) system (4.2) in terms of 
conservation variables. This becomes 



(4.7) 
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( compare (3.15)). Eliminating p t in (4.6) we find that 


J(I + Q) 


P 

pu 

pv 

E 


+ F + G 
x y 


= 0 


_Jt 


(4.8) 


where I is the identity matrix and 


Q = 


0 0 


p 2 2 

R z 2 u -u z 2 -uvz 2 z 2 u 


p 2 2 

R z 2 v -uvz 2 _v z 2 z 2 v 


(4.9) 


R z 3 -uz 3 -vz 3 z 3 . 


and R 2 = (u 2 + v 2 )/2. Comparing (3.14) — (3.17) with (4.6)— (4.9) we see that 
the (p,u,v,p) system leads to a simpler preconditioner than does the 
(p,u,v,S) system. Choi and Merkle [4] pointed out that in the special case 
a = 0, that only the energy equation is modified, i.e., Z£ = 0 when a = 0. 
As before 

-1 8 2 
(I + Q) = I - JL. Q. 

C 


(4.10) 
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5. COMPRESSIBLE ISOENERGETIC SYSTEM 

In the two previous sections we have considered two possibilities for 
adding an additional differential equation to the incompressible equations. A 
different possibility is to use the fact that for the steady state Euler 
equations, when the flow originates from a common reservoir, the specific 
total enthaphy, h = (E + p)/p, is constant throughout the flow. Since we are 
only interested in steady state solutions, we can assume that h = Iiq for all 
time. Such equations have been analyzed by Gottlieb and Gustaffson [7] and 
also Briley, McDonald, and Shararoth [2], Viviand [20], and Rizzi and Eriksson 
[12]. Taking as our unknowns (p,u,v) the equations become in a general 
coordinate system (x,y), 


where 


and 


P*. + -It (Up v + Vp ) + U Y - v X - u Y + v X = 0 
pc 2 pc 2 x y x y x y y x y x 


Ju. + Uu + Vv + (p Y - p Y )/p = 0 (5.1) 

t x y r x y r y x 


Jv + Uv + Vv + (-p X + p x )/p = 0 
t x y *x y *y x 


U = uY - vX , V = -uY + vX , J = X Y - X Y , 
y y x x’ x y y x’ 


P = 


Y ~ 1 


2 . 2 * 
, u + v 

h 0 " 2 


(5. 2 > 


Using a preconditioning similar to that of the previous sections we consider 
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We see that the form of (5.3) is identical to that of (3.4). The only 
difference is that the coefficient p in (3.4) satisfies p = p(p,S) while 
in (5.3) p is given by (5.2). However, the eigenvalue properties of the two 
systems are identical. In particular, it is evident from (5.3) that in the 
absence of preconditioning, i.e., 3 = c and a = 0, that (5.3) is 
simultaneously symmetrizable. In [7] there was an algebraic error and it was 
claimed that the isoenergetic system could not be symmetrized. In [19] we 
presented the matrices that would symmetrize the isoenergetic equations 
written In conservation variables. The proof of symmetry is more obvious 
when ( p , u , v) variables are used as in (5.3). Furthermore, it follows from 
our previous results that (5.3) is symmetrizable for all a and 3 subject 
to the restraint (3.22) 

We can rewrite (5.3) in terras of the conservation variables (p,pu,pv) as 
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z, p + p 
1 *t *t 


J I UZ 2 P t + (pu) t | + F x + G y = °» 


(5.4) 


where 


i vz 2 p t + (pv) . 


[c 2 - + (u 2 + v 2 , 

Bp > ‘ 


and 


. z + P_ q (Y - 1) u 2 + v 2 >, ,2 

2 1 g 2 Y ( h o 2 ) = Z 1 + a/e 


z* “ z 


6 p 


(5.5) 


It follows from (5.2) that 


p = 1 ~ 1 r fu , u 2 + v 2 
Y 


[( h 0 + f~)p fc “ u ^P u ) t “ v(pv) ]. (5.6) 


We can, therefore, rewrite (5.4) so that only time derivatives of 
appear. We thus rewrite (5.4) as 


p, pu, pv 


J(I + Q) ( pu I + F + G 

x y 

,PV/ fc 


where I is the identity matrix and 


= 0 


R 2 z 


1 


Q ' | uR 2 z 2 

2 

WIT z. 


“UZj -vz 


2 

-u z„ 


“UVZ, 


1 


-uvz r 


2 

-v z. 


(5.7) 
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2. U "H V 

where R = ^ also follows that 

-1 B 2 

(I + Q) = I - JL- Q . (5.8) 

c 

As pointed out at the end of Section 3 we can, from (5.6), use p t instead of 
p . Substituting into (5.4) we get a system which is similar to the 
incompressible equations. Now there is no difficulty near stagnation points. 


6. COMPUTATIONAL RESULTS 

Several authors have previously demonstrated the effectiveness of 
preconditioning for the compressible equations with a = 0. Briley, McDonald, 
and Shamroth [2] consider the isoenergetic equations (Section 5) and present 
results using an implicit method for the Navier-Stokes equations with an 
algebraic turbulence model. Rizzi and Eriksson [12] also use a preconditioned 
model based on the analysis of Viviand, also for the isoenergetic equations. 
Rizzi and Eriksson show computer results for the isoenergetic Euler equations 
in both two and three dimensions. These results are calculated using an 
explicit three step Runge-Kutta type algorithm. In addition Choi and Merkle 
[4] analyze a (p,u,v,p) formulation (Section 4) with a = 0, and present 
results for nozzle flow using an implicit A.D.I. type algorithm. In [5] they 
present an alternative approach to preconditioning the equations that is 
effective for very low Mach numbers. 

We present results based on the Runge-Kutta code algorithm described in 
[8], [19]. This is a pseudo-time method that solves the compressible Euler 
equations by a four stage Runge-Kutta formula. The time accuracy is destroyed 
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since a local time step is used together with enthalpy damping. Residual 
smoothing is also applied but without multigrid. We compute flow about a NACA 
0012 with an inflow Mach number at infinity of 0.01 and 0 degree angle of 
attach. Using the original code [8] the residual is reduced two orders of 
magnitued after 1000 Runge-Kutta steps. Using the preconditioning of Section 
3 the residual is reduced by nine orders of magnitude in 1000 steps. In fact, 
the rate of convergence is the same as for transonic cases. 

We note that the analysis presented is based on inviscid flow. Hence, we 
expect the results to also be valid for external high Reynolds number flows. 
However, for flows where the viscous effects are important everywhere, e.g., 
low Reynolds number flows or internal flows, it is not clear that f3 2 should 
vary depending on u 2 + v 2 and not including viscous effects. 
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